home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 6 / Night Owl's Shareware - PDSI-006 - Night Owl Corp (1990).iso / 006a / acm526.zip / ACM526.FOR
Text File  |  1991-02-15  |  61KB  |  1,799 lines
  1. C     ALGORITHM 526, COLLECTED ALGORITHMS FROM ACM. THIS WORK
  2. C     PUBLISHED IN TRANS. MATH. SOFTWARE, 4(2), PP. 160-164.
  3.  
  4. C     PROGRAM  TTIDBS(OUTPUT,TAPE6=OUTPUT)
  5. C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO-
  6. C GRAM PACKAGE.  ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS
  7. C ARE EXPECTED TO BE ZERO.
  8. C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS
  9. C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
  10. C THEREFORE, SYSTEM DEPENDENT.
  11. C DECLARATION STATEMENTS
  12.       DIMENSION   XD(30),YD(30),ZD(30),
  13.      1            XI(6),YI(5),ZI(6,5),
  14.      2            ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5),
  15.      3            IWK(1030),WK(240)
  16.       DATA  NCP/4/
  17.       DATA  NDP/30/
  18.       DATA  XD(1), XD(2), XD(3), XD(4), XD(5), XD(6),
  19.      1      XD(7), XD(8), XD(9), XD(10),XD(11),XD(12),
  20.      2      XD(13),XD(14),XD(15),XD(16),XD(17),XD(18),
  21.      3      XD(19),XD(20),XD(21),XD(22),XD(23),XD(24),
  22.      4      XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/
  23.      5      11.16, 24.20, 19.85, 10.35, 19.72,  0.00,
  24.      6      20.87, 19.99, 10.28,  4.51,  0.00, 16.70,
  25.      7       6.08, 25.00, 14.90,  0.00,  9.66,  5.22,
  26.      8      11.77, 15.10, 25.00, 25.00, 14.59, 15.20,
  27.      9       5.23,  2.14,  0.51, 25.00, 21.67,  3.31/
  28.       DATA  YD(1), YD(2), YD(3), YD(4), YD(5), YD(6),
  29.      1      YD(7), YD(8), YD(9), YD(10),YD(11),YD(12),
  30.      2      YD(13),YD(14),YD(15),YD(16),YD(17),YD(18),
  31.      3      YD(19),YD(20),YD(21),YD(22),YD(23),YD(24),
  32.      4      YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/
  33.      5       1.24, 16.23, 10.72,  4.11,  1.39, 20.00,
  34.      6      20.00,  4.62, 15.16, 20.00,  4.48, 19.65,
  35.      7       4.58, 11.87,  3.12,  0.00, 20.00, 14.66,
  36.      8      10.47, 17.19,  3.87,  0.00,  8.71,  0.00,
  37.      9      10.72, 15.03,  8.37, 20.00, 14.36,  0.13/
  38.       DATA  ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6),
  39.      1      ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12),
  40.      2      ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18),
  41.      3      ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24),
  42.      4      ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/
  43.      5      22.15,  2.83,  7.97, 22.33, 16.83, 34.60,
  44.      6       5.74, 14.72, 21.59, 15.61, 61.77,  6.31,
  45.      7      35.74,  4.40, 21.70, 58.20,  4.73, 40.36,
  46.      8      13.62, 12.57,  8.74, 12.00, 14.81, 21.60,
  47.      9      26.50, 53.10, 49.43,  0.60,  5.52, 44.08/
  48.       DATA  NXI/6/, NYI/5/
  49.       DATA  XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/
  50.      1       0.00,  5.00, 10.00, 15.00, 20.00, 25.00/
  51.       DATA  YI(1), YI(2), YI(3), YI(4), YI(5)/
  52.      1       0.00,  5.00, 10.00, 15.00, 20.00/
  53.       DATA  ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1),
  54.      1      ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2),
  55.      2      ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3),
  56.      3      ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4),
  57.      4      ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/
  58.      5      58.20, 39.55, 26.90, 21.71, 17.68, 12.00,
  59.      6      61.58, 39.39, 22.04, 21.29, 14.36,  8.04,
  60.      7      59.18, 27.39, 16.78, 13.25,  8.59,  5.36,
  61.      8      52.82, 40.27, 22.76, 16.61,  7.40,  2.88,
  62.      9      34.60, 14.05,  4.12,  3.17,  6.31,  0.60/
  63.       DATA  LUN/6/
  64. C CALCULATION
  65.    10 MD=1
  66.       DO 12  IYI=1,NYI
  67.         DO 11  IXI=1,NXI
  68.           IF(IXI.NE.1.OR.IYI.NE.1)  MD=2
  69.           CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI),
  70.      1                ZI1(IXI,IYI),IWK,WK)
  71.    11   CONTINUE
  72.    12 CONTINUE
  73.    15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK)
  74.       DO 17  IYI=1,NYI
  75.         DO 16  IXI=1,NXI
  76.           DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI))
  77.           DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI))
  78.    16   CONTINUE
  79.    17 CONTINUE
  80. C PRINTING OF INPUT DATA
  81.    20 WRITE (LUN,2020)  NDP
  82.       DO 23  IDP=1,NDP
  83.         IF(MOD(IDP,5).EQ.1)    WRITE (LUN,2021)
  84.         WRITE (LUN,2022)  IDP,XD(IDP),YD(IDP),ZD(IDP)
  85.    23 CONTINUE
  86. C PRINTING OF OUTPUT RESULTS
  87.    30 WRITE (LUN,2030)
  88.       WRITE (LUN,2031)  YI
  89.       DO 33  IXI=1,NXI
  90.         WRITE (LUN,2032)  XI(IXI),(ZI1(IXI,IYI),IYI=1,NYI)
  91.    33 CONTINUE
  92.    40 WRITE (LUN,2040)
  93.       WRITE (LUN,2031)  YI
  94.       DO 43  IXI=1,NXI
  95.         WRITE (LUN,2032)  XI(IXI),(DZI1(IXI,IYI),IYI=1,NYI)
  96.    43 CONTINUE
  97.    50 WRITE (LUN,2050)
  98.       WRITE (LUN,2031)  YI
  99.       DO 53  IXI=1,NXI
  100.         WRITE (LUN,2032)  XI(IXI),(ZI2(IXI,IYI),IYI=1,NYI)
  101.    53 CONTINUE
  102.    60 WRITE (LUN,2060)
  103.       WRITE (LUN,2031)  YI
  104.       DO 63  IXI=1,NXI
  105.         WRITE (LUN,2032)  XI(IXI),(DZI2(IXI,IYI),IYI=1,NYI)
  106.    63 CONTINUE
  107.       STOP
  108. C FORMAT STATEMENTS
  109.  2020 FORMAT(1H1,6HTTIDBS/////3X,10HINPUT DATA,8X,5HNDP =,I3///
  110.      1   30H      I      XD     YD     ZD /)
  111.  2021 FORMAT(1X)
  112.  2022 FORMAT(5X,I2,2X,3F7.2)
  113.  2030 FORMAT(1H1,6HTTIDBS/////3X,17HIDBVIP SUBROUTINE///
  114.      1   26X,10HZI1(XI,YI))
  115.  2031 FORMAT(7X,2HXI,4X,3HYI=/12X,5F7.2/)
  116.  2032 FORMAT(1X/1X,F9.2,2X,5F7.2)
  117.  2040 FORMAT(1X/////3X,10HDIFFERENCE///
  118.      1   25X,11HDZI1(XI,YI))
  119.  2050 FORMAT(1H1,6HTTIDBS/////3X,17HIDSFFT SUBROUTINE///
  120.      1   26X,10HZI2(XI,YI))
  121.  2060 FORMAT(1X/////3X,10HDIFFERENCE///
  122.      1   25X,11HDZI2(XI,YI))
  123.       END
  124.       SUBROUTINE  IDBVIP(MD,NCP,NDP,XD,YD,ZD,NIP,XI,YI,ZI,
  125.      1                   IWK,WK)
  126. C THIS SUBROUTINE PERFORMS BIVARIATE INTERPOLATION WHEN THE PRO-
  127. C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
  128. C DISTRIBUTED IN THE PLANE.
  129. C THE INPUT PARAMETERS ARE
  130. C     MD  = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
  131. C         = 1 FOR NEW NCP AND/OR NEW XD-YD,
  132. C         = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
  133. C         = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
  134. C     NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
  135. C           MATING PARTIAL DERIVATIVES AT EACH DATA POINT
  136. C           (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
  137. C     NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
  138. C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE X
  139. C           COORDINATES OF THE DATA POINTS,
  140. C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE Y
  141. C           COORDINATES OF THE DATA POINTS,
  142. C     ZD  = ARRAY OF DIMENSION NDP CONTAINING THE Z
  143. C           COORDINATES OF THE DATA POINTS,
  144. C     NIP = NUMBER OF OUTPUT POINTS AT WHICH INTERPOLATION
  145. C           IS TO BE PERFORMED (MUST BE 1 OR GREATER),
  146. C     XI  = ARRAY OF DIMENSION NIP CONTAINING THE X
  147. C           COORDINATES OF THE OUTPUT POINTS,
  148. C     YI  = ARRAY OF DIMENSION NIP CONTAINING THE Y
  149. C           COORDINATES OF THE OUTPUT POINTS.
  150. C THE OUTPUT PARAMETER IS
  151. C     ZI  = ARRAY OF DIMENSION NIP WHERE INTERPOLATED Z
  152. C           VALUES ARE TO BE STORED.
  153. C THE OTHER PARAMETERS ARE
  154. C     IWK = INTEGER ARRAY OF DIMENSION
  155. C              MAX0(31,27+NCP)*NDP+NIP
  156. C           USED INTERNALLY AS A WORK AREA,
  157. C     WK  = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A
  158. C           WORK AREA.
  159. C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
  160. C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
  161. C YD ARRAYS MUST BE MADE WITH MD=1.  THE CALL WITH MD=2 MUST BE
  162. C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
  163. C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS.  THE CALL WITH
  164. C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
  165. C AND NIP VALUES AND WITH THE SAME CONTENTS OF THE XD, YD, XI,
  166. C AND YI ARRAYS.  BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
  167. C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
  168. C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
  169. C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
  170. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
  171. C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
  172. C THEREFORE, SYSTEM DEPENDENT.
  173. C THIS SUBROUTINE CALLS THE IDCLDP, IDLCTN, IDPDRV, IDPTIP, AND
  174. C IDTANG SUBROUTINES.
  175. C DECLARATION STATEMENTS
  176.       DIMENSION   XD(100),YD(100),ZD(100),XI(1000),YI(1000),
  177.      1            ZI(1000),IWK(4100),WK(800)
  178.       COMMON/IDLC/NIT
  179.       COMMON/IDPI/ITPV
  180.       DATA  LUN/6/
  181. C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
  182. C (FOR MD=1,2,3)
  183.    10 MD0=MD
  184.       NCP0=NCP
  185.       NDP0=NDP
  186.       NIP0=NIP
  187. C ERROR CHECK.  (FOR MD=1,2,3)
  188.    20 IF(MD0.LT.1.OR.MD0.GT.3)           GO TO 90
  189.       IF(NCP0.LT.2.OR.NCP0.GE.NDP0)      GO TO 90
  190.       IF(NDP0.LT.4)                      GO TO 90
  191.       IF(NIP0.LT.1)                      GO TO 90
  192.       IF(MD0.GE.2)        GO TO 21
  193.       IWK(1)=NCP0
  194.       IWK(2)=NDP0
  195.       GO TO 22
  196.    21 NCPPV=IWK(1)
  197.       NDPPV=IWK(2)
  198.       IF(NCP0.NE.NCPPV)   GO TO 90
  199.       IF(NDP0.NE.NDPPV)   GO TO 90
  200.    22 IF(MD0.GE.3)        GO TO 23
  201.       IWK(3)=NIP
  202.       GO TO 30
  203.    23 NIPPV=IWK(3)
  204.       IF(NIP0.NE.NIPPV)   GO TO 90
  205. C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY.  (FOR MD=1,2,3)
  206.    30 JWIPT=16
  207.       JWIWL=6*NDP0+1
  208.       JWIWK=JWIWL
  209.       JWIPL=24*NDP0+1
  210.       JWIWP=30*NDP0+1
  211.       JWIPC=27*NDP0+1
  212.       JWIT0=MAX0(31,27+NCP0)*NDP0
  213. C TRIANGULATES THE X-Y PLANE.  (FOR MD=1)
  214.    40 IF(MD0.GT.1)   GO TO 50
  215.       CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
  216.      1            IWK(JWIWL),IWK(JWIWP),WK)
  217.       IWK(5)=NT
  218.       IWK(6)=NL
  219.       IF(NT.EQ.0)    RETURN
  220. C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT.  (FOR MD=1)
  221.    50 IF(MD0.GT.1)   GO TO 60
  222.       CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC))
  223.       IF(IWK(JWIPC).EQ.0)      RETURN
  224. C LOCATES ALL POINTS AT WHICH INTERPOLATION IS TO BE PERFORMED.
  225. C (FOR MD=1,2)
  226.    60 IF(MD0.EQ.3)   GO TO 70
  227.       NIT=0
  228.       JWIT=JWIT0
  229.       DO 61  IIP=1,NIP0
  230.         JWIT=JWIT+1
  231.         CALL IDLCTN(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
  232.      1            XI(IIP),YI(IIP),IWK(JWIT),IWK(JWIWK),WK)
  233.    61 CONTINUE
  234. C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
  235. C (FOR MD=1,2,3)
  236.    70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
  237. C INTERPOLATES THE ZI VALUES.  (FOR MD=1,2,3)
  238.    80 ITPV=0
  239.       JWIT=JWIT0
  240.       DO 81  IIP=1,NIP0
  241.         JWIT=JWIT+1
  242.         CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
  243.      1              IWK(JWIT),XI(IIP),YI(IIP),ZI(IIP))
  244.    81 CONTINUE
  245.       RETURN
  246. C ERROR EXIT
  247.    90 WRITE (LUN,2090) MD0,NCP0,NDP0,NIP0
  248.       RETURN
  249. C FORMAT STATEMENT FOR ERROR MESSAGE
  250.  2090 FORMAT(1X/41H ***   IMPROPER INPUT PARAMETER VALUE(S)./
  251.      1   7H   MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6,
  252.      2   10X,5HNIP =,I6/
  253.      3   35H ERROR DETECTED IN ROUTINE   IDBVIP/)
  254.       END
  255.       SUBROUTINE  IDCLDP(NDP,XD,YD,NCP,IPC)
  256. C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST
  257. C TO EACH OF THE DATA POINT.
  258. C THE INPUT PARAMETERS ARE
  259. C     NDP = NUMBER OF DATA POINTS,
  260. C     XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
  261. C           COORDINATES OF THE DATA POINTS,
  262. C     NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA
  263. C           POINTS.
  264. C THE OUTPUT PARAMETER IS
  265. C     IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE
  266. C           POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
  267. C           EACH OF THE NDP DATA POINTS ARE TO BE STORED.
  268. C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST
  269. C NOT EXCEED 25.
  270. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
  271. C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
  272. C THEREFORE, SYSTEM DEPENDENT.
  273. C DECLARATION STATEMENTS
  274.       DIMENSION   XD(100),YD(100),IPC(400)
  275.       DIMENSION   DSQ0(25),IPC0(25)
  276.       DATA  NCPMX/25/, LUN/6/
  277. C STATEMENT FUNCTION
  278.       DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
  279. C PRELIMINARY PROCESSING
  280.    10 NDP0=NDP
  281.       NCP0=NCP
  282.       IF(NDP0.LT.2)  GO TO 90
  283.       IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0)    GO TO 90
  284. C CALCULATION
  285.    20 DO 59  IP1=1,NDP0
  286. C - SELECTS NCP POINTS.
  287.         X1=XD(IP1)
  288.         Y1=YD(IP1)
  289.         J1=0
  290.         DSQMX=0.0
  291.         DO 22  IP2=1,NDP0
  292.           IF(IP2.EQ.IP1)  GO TO 22
  293.           DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
  294.           J1=J1+1
  295.           DSQ0(J1)=DSQI
  296.           IPC0(J1)=IP2
  297.           IF(DSQI.LE.DSQMX)    GO TO 21
  298.           DSQMX=DSQI
  299.           JMX=J1
  300.    21     IF(J1.GE.NCP0)  GO TO 23
  301.    22   CONTINUE
  302.    23   IP2MN=IP2+1
  303.         IF(IP2MN.GT.NDP0)      GO TO 30
  304.         DO 25  IP2=IP2MN,NDP0
  305.           IF(IP2.EQ.IP1)  GO TO 25
  306.           DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
  307.           IF(DSQI.GE.DSQMX)    GO TO 25
  308.           DSQ0(JMX)=DSQI
  309.           IPC0(JMX)=IP2
  310.           DSQMX=0.0
  311.           DO 24  J1=1,NCP0
  312.             IF(DSQ0(J1).LE.DSQMX)   GO TO 24
  313.             DSQMX=DSQ0(J1)
  314.             JMX=J1
  315.    24     CONTINUE
  316.    25   CONTINUE
  317. C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR.
  318.    30   IP2=IPC0(1)
  319.         DX12=XD(IP2)-X1
  320.         DY12=YD(IP2)-Y1
  321.         DO 31  J3=2,NCP0
  322.           IP3=IPC0(J3)
  323.           DX13=XD(IP3)-X1
  324.           DY13=YD(IP3)-Y1
  325.           IF((DY13*DX12-DX13*DY12).NE.0.0)    GO TO 50
  326.    31   CONTINUE
  327. C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT.
  328.    40   NCLPT=0
  329.         DO 43  IP3=1,NDP0
  330.           IF(IP3.EQ.IP1)       GO TO 43
  331.           DO 41  J4=1,NCP0
  332.             IF(IP3.EQ.IPC0(J4))     GO TO 43
  333.    41     CONTINUE
  334.           DX13=XD(IP3)-X1
  335.           DY13=YD(IP3)-Y1
  336.           IF((DY13*DX12-DX13*DY12).EQ.0.0)    GO TO 43
  337.           DSQI=DSQF(X1,Y1,XD(IP3),YD(IP3))
  338.           IF(NCLPT.EQ.0)       GO TO 42
  339.           IF(DSQI.GE.DSQMN)    GO TO 43
  340.    42     NCLPT=1
  341.           DSQMN=DSQI
  342.           IP3MN=IP3
  343.    43   CONTINUE
  344.         IF(NCLPT.EQ.0)    GO TO 91
  345.         DSQMX=DSQMN
  346.         IPC0(JMX)=IP3MN
  347. C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY.
  348.    50   J1=(IP1-1)*NCP0
  349.         DO 51  J2=1,NCP0
  350.           J1=J1+1
  351.           IPC(J1)=IPC0(J2)
  352.    51   CONTINUE
  353.    59 CONTINUE
  354.       RETURN
  355. C ERROR EXIT
  356.    90 WRITE (LUN,2090)
  357.       GO TO 92
  358.    91 WRITE (LUN,2091)
  359.    92 WRITE (LUN,2092)  NDP0,NCP0
  360.       IPC(1)=0
  361.       RETURN
  362. C FORMAT STATEMENTS FOR ERROR MESSAGES
  363.  2090 FORMAT(1X/41H ***   IMPROPER INPUT PARAMETER VALUE(S).)
  364.  2091 FORMAT(1X/33H ***   ALL COLLINEAR DATA POINTS.)
  365.  2092 FORMAT(8H   NDP =,I5,5X,5HNCP =,I5/
  366.      1   35H ERROR DETECTED IN ROUTINE   IDCLDP/)
  367.       END
  368.       SUBROUTINE IDGRID(XD, YD, NT, IPT, NL, IPL, NXI, NYI, XI, YI,
  369.      *  NGP, IGP)
  370. C THIS SUBROUTINE ORGANIZES GRID POINTS FOR SURFACE FITTING BY
  371. C SORTING THEM IN ASCENDING ORDER OF TRIANGLE NUMBERS AND OF THE
  372. C BORDER LINE SEGMENT NUMBER.
  373. C THE INPUT PARAMETERS ARE
  374. C     XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
  375. C           COORDINATES OF THE DATA POINTS, WHERE NDP IS THE
  376. C           NUMBER OF THE DATA POINTS,
  377. C     NT  = NUMBER OF TRIANGLES,
  378. C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
  379. C           POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
  380. C     NL  = NUMBER OF BORDER LINE SEGMENTS,
  381. C     IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
  382. C           POINT NUMBERS OF THE END POINTS OF THE BORDER
  383. C           LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
  384. C           NUMBERS,
  385. C     NXI = NUMBER OF GRID POINTS IN THE X COORDINATE,
  386. C     NYI = NUMBER OF GRID POINTS IN THE Y COORDINATE,
  387. C     XI,YI = ARRAYS OF DIMENSION NXI AND NYI CONTAINING
  388. C           THE X AND Y COORDINATES OF THE GRID POINTS,
  389. C           RESPECTIVELY.
  390. C THE OUTPUT PARAMETERS ARE
  391. C     NGP = INTEGER ARRAY OF DIMENSION 2*(NT+2*NL) WHERE THE
  392. C           NUMBER OF GRID POINTS THAT BELONG TO EACH OF THE
  393. C           TRIANGLES OR OF THE BORDER LINE SEGMENTS ARE TO
  394. C           BE STORED,
  395. C     IGP = INTEGER ARRAY OF DIMENSION NXI*NYI WHERE THE
  396. C           GRID POINT NUMBERS ARE TO BE STORED IN ASCENDING
  397. C           ORDER OF THE TRIANGLE NUMBER AND THE BORDER LINE
  398. C           SEGMENT NUMBER.
  399. C DECLARATION STATEMENTS
  400.       DIMENSION XD(100), YD(100), IPT(585), IPL(300), XI(101),
  401.      *  YI(101), NGP(800), IGP(10201)
  402. C STATEMENT FUNCTIONS
  403.       SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3)
  404.       SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2)
  405. C PRELIMINARY PROCESSING
  406.       NT0 = NT
  407.       NL0 = NL
  408.       NXI0 = NXI
  409.       NYI0 = NYI
  410.       NXINYI = NXI0*NYI0
  411.       XIMN = AMIN1(XI(1),XI(NXI0))
  412.       XIMX = AMAX1(XI(1),XI(NXI0))
  413.       YIMN = AMIN1(YI(1),YI(NYI0))
  414.       YIMX = AMAX1(YI(1),YI(NYI0))
  415. C DETERMINES GRID POINTS INSIDE THE DATA AREA.
  416.       JNGP0 = 0
  417.       JNGP1 = 2*(NT0+2*NL0) + 1
  418.       JIGP0 = 0
  419.       JIGP1 = NXINYI + 1
  420.       DO 160 IT0=1,NT0
  421.         NGP0 = 0
  422.         NGP1 = 0
  423.         IT0T3 = IT0*3
  424.         IP1 = IPT(IT0T3-2)
  425.         IP2 = IPT(IT0T3-1)
  426.         IP3 = IPT(IT0T3)
  427.         X1 = XD(IP1)
  428.         Y1 = YD(IP1)
  429.         X2 = XD(IP2)
  430.         Y2 = YD(IP2)
  431.         X3 = XD(IP3)
  432.         Y3 = YD(IP3)
  433.         XMN = AMIN1(X1,X2,X3)
  434.         XMX = AMAX1(X1,X2,X3)
  435.         YMN = AMIN1(Y1,Y2,Y3)
  436.         YMX = AMAX1(Y1,Y2,Y3)
  437.         INSD = 0
  438.         DO 20 IXI=1,NXI0
  439.           IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 10
  440.           IF (INSD.EQ.0) GO TO 20
  441.           IXIMX = IXI - 1
  442.           GO TO 30
  443.    10     IF (INSD.EQ.1) GO TO 20
  444.           INSD = 1
  445.           IXIMN = IXI
  446.    20   CONTINUE
  447.         IF (INSD.EQ.0) GO TO 150
  448.         IXIMX = NXI0
  449.    30   DO 140 IYI=1,NYI0
  450.           YII = YI(IYI)
  451.           IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 140
  452.           DO 130 IXI=IXIMN,IXIMX
  453.             XII = XI(IXI)
  454.             L = 0
  455.             IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 130, 40, 50
  456.    40       L = 1
  457.    50       IF (SIDE(X2,Y2,X3,Y3,XII,YII)) 130, 60, 70
  458.    60       L = 1
  459.    70       IF (SIDE(X3,Y3,X1,Y1,XII,YII)) 130, 80, 90
  460.    80       L = 1
  461.    90       IZI = NXI0*(IYI-1) + IXI
  462.             IF (L.EQ.1) GO TO 100
  463.             NGP0 = NGP0 + 1
  464.             JIGP0 = JIGP0 + 1
  465.             IGP(JIGP0) = IZI
  466.             GO TO 130
  467.   100       IF (JIGP1.GT.NXINYI) GO TO 120
  468.             DO 110 JIGP1I=JIGP1,NXINYI
  469.               IF (IZI.EQ.IGP(JIGP1I)) GO TO 130
  470.   110       CONTINUE
  471.   120       NGP1 = NGP1 + 1
  472.             JIGP1 = JIGP1 - 1
  473.             IGP(JIGP1) = IZI
  474.   130     CONTINUE
  475.   140   CONTINUE
  476.   150   JNGP0 = JNGP0 + 1
  477.         NGP(JNGP0) = NGP0
  478.         JNGP1 = JNGP1 - 1
  479.         NGP(JNGP1) = NGP1
  480.   160 CONTINUE
  481. C DETERMINES GRID POINTS OUTSIDE THE DATA AREA.
  482. C - IN SEMI-INFINITE RECTANGULAR AREA.
  483.       DO 450 IL0=1,NL0
  484.         NGP0 = 0
  485.         NGP1 = 0
  486.         IL0T3 = IL0*3
  487.         IP1 = IPL(IL0T3-2)
  488.         IP2 = IPL(IL0T3-1)
  489.         X1 = XD(IP1)
  490.         Y1 = YD(IP1)
  491.         X2 = XD(IP2)
  492.         Y2 = YD(IP2)
  493.         XMN = XIMN
  494.         XMX = XIMX
  495.         YMN = YIMN
  496.         YMX = YIMX
  497.         IF (Y2.GE.Y1) XMN = AMIN1(X1,X2)
  498.         IF (Y2.LE.Y1) XMX = AMAX1(X1,X2)
  499.         IF (X2.LE.X1) YMN = AMIN1(Y1,Y2)
  500.         IF (X2.GE.X1) YMX = AMAX1(Y1,Y2)
  501.         INSD = 0
  502.         DO 180 IXI=1,NXI0
  503.           IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 170
  504.           IF (INSD.EQ.0) GO TO 180
  505.           IXIMX = IXI - 1
  506.           GO TO 190
  507.   170     IF (INSD.EQ.1) GO TO 180
  508.           INSD = 1
  509.           IXIMN = IXI
  510.   180   CONTINUE
  511.         IF (INSD.EQ.0) GO TO 310
  512.         IXIMX = NXI0
  513.   190   DO 300 IYI=1,NYI0
  514.           YII = YI(IYI)
  515.           IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 300
  516.           DO 290 IXI=IXIMN,IXIMX
  517.             XII = XI(IXI)
  518.             L = 0
  519.             IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 210, 200, 290
  520.   200       L = 1
  521.   210       IF (SPDT(X2,Y2,X1,Y1,XII,YII)) 290, 220, 230
  522.   220       L = 1
  523.   230       IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 290, 240, 250
  524.   240       L = 1
  525.   250       IZI = NXI0*(IYI-1) + IXI
  526.             IF (L.EQ.1) GO TO 260
  527.             NGP0 = NGP0 + 1
  528.             JIGP0 = JIGP0 + 1
  529.             IGP(JIGP0) = IZI
  530.             GO TO 290
  531.   260       IF (JIGP1.GT.NXINYI) GO TO 280
  532.             DO 270 JIGP1I=JIGP1,NXINYI
  533.               IF (IZI.EQ.IGP(JIGP1I)) GO TO 290
  534.   270       CONTINUE
  535.   280       NGP1 = NGP1 + 1
  536.             JIGP1 = JIGP1 - 1
  537.             IGP(JIGP1) = IZI
  538.   290     CONTINUE
  539.   300   CONTINUE
  540.   310   JNGP0 = JNGP0 + 1
  541.         NGP(JNGP0) = NGP0
  542.         JNGP1 = JNGP1 - 1
  543.         NGP(JNGP1) = NGP1
  544. C - IN SEMI-INFINITE TRIANGULAR AREA.
  545.         NGP0 = 0
  546.         NGP1 = 0
  547.         ILP1 = MOD(IL0,NL0) + 1
  548.         ILP1T3 = ILP1*3
  549.         IP3 = IPL(ILP1T3-1)
  550.         X3 = XD(IP3)
  551.         Y3 = YD(IP3)
  552.         XMN = XIMN
  553.         XMX = XIMX
  554.         YMN = YIMN
  555.         YMX = YIMX
  556.         IF (Y3.GE.Y2 .AND. Y2.GE.Y1) XMN = X2
  557.         IF (Y3.LE.Y2 .AND. Y2.LE.Y1) XMX = X2
  558.         IF (X3.LE.X2 .AND. X2.LE.X1) YMN = Y2
  559.         IF (X3.GE.X2 .AND. X2.GE.X1) YMX = Y2
  560.         INSD = 0
  561.         DO 330 IXI=1,NXI0
  562.           IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 320
  563.           IF (INSD.EQ.0) GO TO 330
  564.           IXIMX = IXI - 1
  565.           GO TO 340
  566.   320     IF (INSD.EQ.1) GO TO 330
  567.           INSD = 1
  568.           IXIMN = IXI
  569.   330   CONTINUE
  570.         IF (INSD.EQ.0) GO TO 440
  571.         IXIMX = NXI0
  572.   340   DO 430 IYI=1,NYI0
  573.           YII = YI(IYI)
  574.           IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 430
  575.           DO 420 IXI=IXIMN,IXIMX
  576.             XII = XI(IXI)
  577.             L = 0
  578.             IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 360, 350, 420
  579.   350       L = 1
  580.   360       IF (SPDT(X3,Y3,X2,Y2,XII,YII)) 380, 370, 420
  581.   370       L = 1
  582.   380       IZI = NXI0*(IYI-1) + IXI
  583.             IF (L.EQ.1) GO TO 390
  584.             NGP0 = NGP0 + 1
  585.             JIGP0 = JIGP0 + 1
  586.             IGP(JIGP0) = IZI
  587.             GO TO 420
  588.   390       IF (JIGP1.GT.NXINYI) GO TO 410
  589.             DO 400 JIGP1I=JIGP1,NXINYI
  590.               IF (IZI.EQ.IGP(JIGP1I)) GO TO 420
  591.   400       CONTINUE
  592.   410       NGP1 = NGP1 + 1
  593.             JIGP1 = JIGP1 - 1
  594.             IGP(JIGP1) = IZI
  595.   420     CONTINUE
  596.   430   CONTINUE
  597.   440   JNGP0 = JNGP0 + 1
  598.         NGP(JNGP0) = NGP0
  599.         JNGP1 = JNGP1 - 1
  600.         NGP(JNGP1) = NGP1
  601.   450 CONTINUE
  602.       RETURN
  603.       END
  604.       SUBROUTINE IDLCTN(NDP, XD, YD, NT, IPT, NL, IPL, XII, YII, ITI,
  605.      *  IWK, WK)
  606. C THIS SUBROUTINE LOCATES A POINT, I.E., DETERMINES TO WHAT TRI-
  607. C ANGLE A GIVEN POINT (XII,YII) BELONGS.  WHEN THE GIVEN POINT
  608. C DOES NOT LIE INSIDE THE DATA AREA, THIS SUBROUTINE DETERMINES
  609. C THE BORDER LINE SEGMENT WHEN THE POINT LIES IN AN OUTSIDE
  610. C RECTANGULAR AREA, AND TWO BORDER LINE SEGMENTS WHEN THE POINT
  611. C LIES IN AN OUTSIDE TRIANGULAR AREA.
  612. C THE INPUT PARAMETERS ARE
  613. C     NDP = NUMBER OF DATA POINTS,
  614. C     XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
  615. C           COORDINATES OF THE DATA POINTS,
  616. C     NT  = NUMBER OF TRIANGLES,
  617. C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
  618. C           POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
  619. C     NL  = NUMBER OF BORDER LINE SEGMENTS,
  620. C     IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
  621. C           POINT NUMBERS OF THE END POINTS OF THE BORDER
  622. C           LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
  623. C           NUMBERS,
  624. C     XII,YII = X AND Y COORDINATES OF THE POINT TO BE
  625. C           LOCATED.
  626. C THE OUTPUT PARAMETER IS
  627. C     ITI = TRIANGLE NUMBER, WHEN THE POINT IS INSIDE THE
  628. C           DATA AREA, OR
  629. C           TWO BORDER LINE SEGMENT NUMBERS, IL1 AND IL2,
  630. C           CODED TO IL1*(NT+NL)+IL2, WHEN THE POINT IS
  631. C           OUTSIDE THE DATA AREA.
  632. C THE OTHER PARAMETERS ARE
  633. C     IWK = INTEGER ARRAY OF DIMENSION 18*NDP USED INTER-
  634. C           NALLY AS A WORK AREA,
  635. C     WK  = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A
  636. C           WORK AREA.
  637. C DECLARATION STATEMENTS
  638.       DIMENSION XD(100), YD(100), IPT(585), IPL(300), IWK(1800),
  639.      *  WK(800)
  640.       DIMENSION NTSC(9), IDSC(9)
  641.       COMMON /IDLC/ NIT
  642. C STATEMENT FUNCTIONS
  643.       SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3)
  644.       SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2)
  645. C PRELIMINARY PROCESSING
  646.       NDP0 = NDP
  647.       NT0 = NT
  648.       NL0 = NL
  649.       NTL = NT0 + NL0
  650.       X0 = XII
  651.       Y0 = YII
  652. C PROCESSING FOR A NEW SET OF DATA POINTS
  653.       IF (NIT.NE.0) GO TO 80
  654.       NIT = 1
  655. C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS.
  656.       XMN = XD(1)
  657.       XMX = XMN
  658.       YMN = YD(1)
  659.       YMX = YMN
  660.       DO 10 IDP=2,NDP0
  661.         XI = XD(IDP)
  662.         YI = YD(IDP)
  663.         XMN = AMIN1(XI,XMN)
  664.         XMX = AMAX1(XI,XMX)
  665.         YMN = AMIN1(YI,YMN)
  666.         YMX = AMAX1(YI,YMX)
  667.    10 CONTINUE
  668.       XS1 = (XMN+XMN+XMX)/3.0
  669.       XS2 = (XMN+XMX+XMX)/3.0
  670.       YS1 = (YMN+YMN+YMX)/3.0
  671.       YS2 = (YMN+YMX+YMX)/3.0
  672. C - DETERMINES AND STORES IN THE IWK ARRAY TRIANGLE NUMBERS OF
  673. C - THE TRIANGLES ASSOCIATED WITH EACH OF THE NINE SECTIONS.
  674.       DO 20 ISC=1,9
  675.         NTSC(ISC) = 0
  676.         IDSC(ISC) = 0
  677.    20 CONTINUE
  678.       IT0T3 = 0
  679.       JWK = 0
  680.       DO 70 IT0=1,NT0
  681.         IT0T3 = IT0T3 + 3
  682.         I1 = IPT(IT0T3-2)
  683.         I2 = IPT(IT0T3-1)
  684.         I3 = IPT(IT0T3)
  685.         XMN = AMIN1(XD(I1),XD(I2),XD(I3))
  686.         XMX = AMAX1(XD(I1),XD(I2),XD(I3))
  687.         YMN = AMIN1(YD(I1),YD(I2),YD(I3))
  688.         YMX = AMAX1(YD(I1),YD(I2),YD(I3))
  689.         IF (YMN.GT.YS1) GO TO 30
  690.         IF (XMN.LE.XS1) IDSC(1) = 1
  691.         IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(2) = 1
  692.         IF (XMX.GE.XS2) IDSC(3) = 1
  693.    30   IF (YMX.LT.YS1 .OR. YMN.GT.YS2) GO TO 40
  694.         IF (XMN.LE.XS1) IDSC(4) = 1
  695.         IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(5) = 1
  696.         IF (XMX.GE.XS2) IDSC(6) = 1
  697.    40   IF (YMX.LT.YS2) GO TO 50
  698.         IF (XMN.LE.XS1) IDSC(7) = 1
  699.         IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(8) = 1
  700.         IF (XMX.GE.XS2) IDSC(9) = 1
  701.    50   DO 60 ISC=1,9
  702.           IF (IDSC(ISC).EQ.0) GO TO 60
  703.           JIWK = 9*NTSC(ISC) + ISC
  704.           IWK(JIWK) = IT0
  705.           NTSC(ISC) = NTSC(ISC) + 1
  706.           IDSC(ISC) = 0
  707.    60   CONTINUE
  708. C - STORES IN THE WK ARRAY THE MINIMUM AND MAXIMUM OF THE X AND
  709. C - Y COORDINATE VALUES FOR EACH OF THE TRIANGLE.
  710.         JWK = JWK + 4
  711.         WK(JWK-3) = XMN
  712.         WK(JWK-2) = XMX
  713.         WK(JWK-1) = YMN
  714.         WK(JWK) = YMX
  715.    70 CONTINUE
  716.       GO TO 110
  717. C CHECKS IF IN THE SAME TRIANGLE AS PREVIOUS.
  718.    80 IT0 = ITIPV
  719.       IF (IT0.GT.NT0) GO TO 90
  720.       IT0T3 = IT0*3
  721.       IP1 = IPT(IT0T3-2)
  722.       X1 = XD(IP1)
  723.       Y1 = YD(IP1)
  724.       IP2 = IPT(IT0T3-1)
  725.       X2 = XD(IP2)
  726.       Y2 = YD(IP2)
  727.       IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110
  728.       IP3 = IPT(IT0T3)
  729.       X3 = XD(IP3)
  730.       Y3 = YD(IP3)
  731.       IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 110
  732.       IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 110
  733.       GO TO 170
  734. C CHECKS IF ON THE SAME BORDER LINE SEGMENT.
  735.    90 IL1 = IT0/NTL
  736.       IL2 = IT0 - IL1*NTL
  737.       IL1T3 = IL1*3
  738.       IP1 = IPL(IL1T3-2)
  739.       X1 = XD(IP1)
  740.       Y1 = YD(IP1)
  741.       IP2 = IPL(IL1T3-1)
  742.       X2 = XD(IP2)
  743.       Y2 = YD(IP2)
  744.       IF (IL2.NE.IL1) GO TO 100
  745.       IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110
  746.       IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 110
  747.       IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110
  748.       GO TO 170
  749. C CHECKS IF BETWEEN THE SAME TWO BORDER LINE SEGMENTS.
  750.   100 IF (SPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110
  751.       IP3 = IPL(3*IL2-1)
  752.       X3 = XD(IP3)
  753.       Y3 = YD(IP3)
  754.       IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 170
  755. C LOCATES INSIDE THE DATA AREA.
  756. C - DETERMINES THE SECTION IN WHICH THE POINT IN QUESTION LIES.
  757.   110 ISC = 1
  758.       IF (X0.GE.XS1) ISC = ISC + 1
  759.       IF (X0.GE.XS2) ISC = ISC + 1
  760.       IF (Y0.GE.YS1) ISC = ISC + 3
  761.       IF (Y0.GE.YS2) ISC = ISC + 3
  762. C - SEARCHES THROUGH THE TRIANGLES ASSOCIATED WITH THE SECTION.
  763.       NTSCI = NTSC(ISC)
  764.       IF (NTSCI.LE.0) GO TO 130
  765.       JIWK = -9 + ISC
  766.       DO 120 ITSC=1,NTSCI
  767.         JIWK = JIWK + 9
  768.         IT0 = IWK(JIWK)
  769.         JWK = IT0*4
  770.         IF (X0.LT.WK(JWK-3)) GO TO 120
  771.         IF (X0.GT.WK(JWK-2)) GO TO 120
  772.         IF (Y0.LT.WK(JWK-1)) GO TO 120
  773.         IF (Y0.GT.WK(JWK)) GO TO 120
  774.         IT0T3 = IT0*3
  775.         IP1 = IPT(IT0T3-2)
  776.         X1 = XD(IP1)
  777.         Y1 = YD(IP1)
  778.         IP2 = IPT(IT0T3-1)
  779.         X2 = XD(IP2)
  780.         Y2 = YD(IP2)
  781.         IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 120
  782.         IP3 = IPT(IT0T3)
  783.         X3 = XD(IP3)
  784.         Y3 = YD(IP3)
  785.         IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 120
  786.         IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 120
  787.         GO TO 170
  788.   120 CONTINUE
  789. C LOCATES OUTSIDE THE DATA AREA.
  790.   130 DO 150 IL1=1,NL0
  791.         IL1T3 = IL1*3
  792.         IP1 = IPL(IL1T3-2)
  793.         X1 = XD(IP1)
  794.         Y1 = YD(IP1)
  795.         IP2 = IPL(IL1T3-1)
  796.         X2 = XD(IP2)
  797.         Y2 = YD(IP2)
  798.         IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 150
  799.         IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 140
  800.         IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 150
  801.         IL2 = IL1
  802.         GO TO 160
  803.   140   IL2 = MOD(IL1,NL0) + 1
  804.         IP3 = IPL(3*IL2-1)
  805.         X3 = XD(IP3)
  806.         Y3 = YD(IP3)
  807.         IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 160
  808.   150 CONTINUE
  809.       IT0 = 1
  810.       GO TO 170
  811.   160 IT0 = IL1*NTL + IL2
  812. C NORMAL EXIT
  813.   170 ITI = IT0
  814.       ITIPV = IT0
  815.       RETURN
  816.       END
  817.       SUBROUTINE  IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD)
  818. C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND
  819. C SECOND ORDER AT THE DATA POINTS.
  820. C THE INPUT PARAMETERS ARE
  821. C     NDP = NUMBER OF DATA POINTS,
  822. C     XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
  823. C           Y, AND Z COORDINATES OF THE DATA POINTS,
  824. C     NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
  825. C           MATING PARTIAL DERIVATIVES AT EACH DATA POINT,
  826. C     IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING
  827. C           THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
  828. C           EACH OF THE NDP DATA POINTS.
  829. C THE OUTPUT PARAMETER IS
  830. C     PD  = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED
  831. C           ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA
  832. C           POINTS ARE TO BE STORED.
  833. C DECLARATION STATEMENTS
  834.       DIMENSION   XD(100),YD(100),ZD(100),IPC(400),PD(500)
  835.       REAL        NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY
  836. C PRELIMINARY PROCESSING
  837.    10 NDP0=NDP
  838.       NCP0=NCP
  839.       NCPM1=NCP0-1
  840. C ESTIMATION OF ZX AND ZY
  841.   20  DO 24  IP0=1,NDP0
  842.         X0=XD(IP0)
  843.         Y0=YD(IP0)
  844.         Z0=ZD(IP0)
  845.         NMX=0.0
  846.         NMY=0.0
  847.         NMZ=0.0
  848.         JIPC0=NCP0*(IP0-1)
  849.         DO 23  IC1=1,NCPM1
  850.           JIPC=JIPC0+IC1
  851.           IPI=IPC(JIPC)
  852.           DX1=XD(IPI)-X0
  853.           DY1=YD(IPI)-Y0
  854.           DZ1=ZD(IPI)-Z0
  855.           IC2MN=IC1+1
  856.           DO 22  IC2=IC2MN,NCP0
  857.             JIPC=JIPC0+IC2
  858.             IPI=IPC(JIPC)
  859.             DX2=XD(IPI)-X0
  860.             DY2=YD(IPI)-Y0
  861.             DNMZ=DX1*DY2-DY1*DX2
  862.             IF(DNMZ.EQ.0.0)    GO TO 22
  863.             DZ2=ZD(IPI)-Z0
  864.             DNMX=DY1*DZ2-DZ1*DY2
  865.             DNMY=DZ1*DX2-DX1*DZ2
  866.             IF(DNMZ.GE.0.0)    GO TO 21
  867.             DNMX=-DNMX
  868.             DNMY=-DNMY
  869.             DNMZ=-DNMZ
  870.    21       NMX=NMX+DNMX
  871.             NMY=NMY+DNMY
  872.             NMZ=NMZ+DNMZ
  873.    22     CONTINUE
  874.    23   CONTINUE
  875.         JPD0=5*IP0
  876.         PD(JPD0-4)=-NMX/NMZ
  877.         PD(JPD0-3)=-NMY/NMZ
  878.    24 CONTINUE
  879. C ESTIMATION OF ZXX, ZXY, AND ZYY
  880.    30 DO 34  IP0=1,NDP0
  881.         JPD0=JPD0+5
  882.         X0=XD(IP0)
  883.         JPD0=5*IP0
  884.         Y0=YD(IP0)
  885.         ZX0=PD(JPD0-4)
  886.         ZY0=PD(JPD0-3)
  887.         NMXX=0.0
  888.         NMXY=0.0
  889.         NMYX=0.0
  890.         NMYY=0.0
  891.         NMZ =0.0
  892.         JIPC0=NCP0*(IP0-1)
  893.         DO 33  IC1=1,NCPM1
  894.           JIPC=JIPC0+IC1
  895.           IPI=IPC(JIPC)
  896.           DX1=XD(IPI)-X0
  897.           DY1=YD(IPI)-Y0
  898.           JPD=5*IPI
  899.           DZX1=PD(JPD-4)-ZX0
  900.           DZY1=PD(JPD-3)-ZY0
  901.           IC2MN=IC1+1
  902.           DO 32  IC2=IC2MN,NCP0
  903.             JIPC=JIPC0+IC2
  904.             IPI=IPC(JIPC)
  905.             DX2=XD(IPI)-X0
  906.             DY2=YD(IPI)-Y0
  907.             DNMZ =DX1*DY2 -DY1*DX2
  908.             IF(DNMZ.EQ.0.0)    GO TO 32
  909.             JPD=5*IPI
  910.             DZX2=PD(JPD-4)-ZX0
  911.             DZY2=PD(JPD-3)-ZY0
  912.             DNMXX=DY1*DZX2-DZX1*DY2
  913.             DNMXY=DZX1*DX2-DX1*DZX2
  914.             DNMYX=DY1*DZY2-DZY1*DY2
  915.             DNMYY=DZY1*DX2-DX1*DZY2
  916.             IF(DNMZ.GE.0.0)    GO TO 31
  917.             DNMXX=-DNMXX
  918.             DNMXY=-DNMXY
  919.             DNMYX=-DNMYX
  920.             DNMYY=-DNMYY
  921.             DNMZ =-DNMZ
  922.    31       NMXX=NMXX+DNMXX
  923.             NMXY=NMXY+DNMXY
  924.             NMYX=NMYX+DNMYX                                             
  925. ID010010
  926.             NMYY=NMYY+DNMYY
  927.             NMZ =NMZ +DNMZ
  928.    32     CONTINUE
  929.    33   CONTINUE
  930.         PD(JPD0-2)=-NMXX/NMZ
  931.         PD(JPD0-1)=-(NMXY+NMYX)/(2.0*NMZ)
  932.         PD(JPD0)  =-NMYY/NMZ
  933.    34 CONTINUE
  934.       RETURN
  935.       END
  936.       SUBROUTINE  IDPTIP(XD,YD,ZD,NT,IPT,NL,IPL,PDD,ITI,XII,YII,
  937.      1                   ZII)
  938. C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPOLA-
  939. C TION, I.E., DETERMINES THE Z VALUE AT A POINT.
  940. C THE INPUT PARAMETERS ARE
  941. C     XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
  942. C           Y, AND Z COORDINATES OF THE DATA POINTS, WHERE
  943. C           NDP IS THE NUMBER OF THE DATA POINTS,
  944. C     NT  = NUMBER OF TRIANGLES,
  945. C     IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE
  946. C           POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES,
  947. C     NL  = NUMBER OF BORDER LINE SEGMENTS,
  948. C     IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE
  949. C           POINT NUMBERS OF THE END POINTS OF THE BORDER
  950. C           LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE
  951. C           NUMBERS,
  952. C     PDD = ARRAY OF DIMENSION 5*NDP CONTAINING THE PARTIAL
  953. C           DERIVATIVES AT THE DATA POINTS,
  954. C     ITI = TRIANGLE NUMBER OF THE TRIANGLE IN WHICH LIES
  955. C           THE POINT FOR WHICH INTERPOLATION IS TO BE
  956. C           PERFORMED,
  957. C     XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH
  958. C           INTERPOLATION IS TO BE PERFORMED.
  959. C THE OUTPUT PARAMETER IS
  960. C     ZII = INTERPOLATED Z VALUE.
  961. C DECLARATION STATEMENTS
  962.       DIMENSION   XD(100),YD(100),ZD(100),IPT(585),IPL(300),
  963.      1            PDD(500)
  964.       COMMON/IDPI/ITPV
  965.       DIMENSION   X(3),Y(3),Z(3),PD(15),
  966.      1            ZU(3),ZV(3),ZUU(3),ZUV(3),ZVV(3)
  967.       REAL        LU,LV
  968.       EQUIVALENCE (P5,P50)
  969. C PRELIMINARY PROCESSING
  970.    10 IT0=ITI
  971.       NTL=NT+NL
  972.       IF(IT0.LE.NTL)      GO TO 20
  973.       IL1=IT0/NTL
  974.       IL2=IT0-IL1*NTL
  975.       IF(IL1.EQ.IL2)      GO TO 40
  976.       GO TO 60
  977. C CALCULATION OF ZII BY INTERPOLATION.
  978. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
  979.    20 IF(IT0.EQ.ITPV)     GO TO 30
  980. C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE
  981. C VERTEXES.
  982.    21 JIPT=3*(IT0-1)
  983.       JPD=0
  984.       DO 23  I=1,3
  985.         JIPT=JIPT+1
  986.         IDP=IPT(JIPT)
  987.         X(I)=XD(IDP)
  988.         Y(I)=YD(IDP)
  989.         Z(I)=ZD(IDP)
  990.         JPDD=5*(IDP-1)
  991.         DO 22  KPD=1,5
  992.           JPD=JPD+1
  993.           JPDD=JPDD+1
  994.           PD(JPD)=PDD(JPDD)
  995.    22   CONTINUE
  996.    23 CONTINUE
  997. C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
  998. C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
  999. C AND VICE VERSA.
  1000.    24 X0=X(1)
  1001.       Y0=Y(1)
  1002.       A=X(2)-X0
  1003.       B=X(3)-X0
  1004.       C=Y(2)-Y0
  1005.       D=Y(3)-Y0
  1006.       AD=A*D
  1007.       BC=B*C
  1008.       DLT=AD-BC
  1009.       AP= D/DLT
  1010.       BP=-B/DLT
  1011.       CP=-C/DLT
  1012.       DP= A/DLT
  1013. C CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE
  1014. C TRIANGLE FOR THE U-V COORDINATE SYSTEM.
  1015.    25 AA=A*A
  1016.       ACT2=2.0*A*C
  1017.       CC=C*C
  1018.       AB=A*B
  1019.       ADBC=AD+BC
  1020.       CD=C*D
  1021.       BB=B*B
  1022.       BDT2=2.0*B*D
  1023.       DD=D*D
  1024.       DO 26  I=1,3
  1025.         JPD=5*I
  1026.         ZU(I)=A*PD(JPD-4)+C*PD(JPD-3)
  1027.         ZV(I)=B*PD(JPD-4)+D*PD(JPD-3)
  1028.         ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD)
  1029.         ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD)
  1030.         ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD)
  1031.    26 CONTINUE
  1032. C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
  1033.    27 P00=Z(1)
  1034.       P10=ZU(1)
  1035.       P01=ZV(1)
  1036.       P20=0.5*ZUU(1)
  1037.       P11=ZUV(1)
  1038.       P02=0.5*ZVV(1)
  1039.       H1=Z(2)-P00-P10-P20
  1040.       H2=ZU(2)-P10-ZUU(1)
  1041.       H3=ZUU(2)-ZUU(1)
  1042.       P30= 10.0*H1-4.0*H2+0.5*H3
  1043.       P40=-15.0*H1+7.0*H2    -H3
  1044.       P50=  6.0*H1-3.0*H2+0.5*H3
  1045.       H1=Z(3)-P00-P01-P02
  1046.       H2=ZV(3)-P01-ZVV(1)
  1047.       H3=ZVV(3)-ZVV(1)
  1048.       P03= 10.0*H1-4.0*H2+0.5*H3
  1049.       P04=-15.0*H1+7.0*H2    -H3
  1050.       P05=  6.0*H1-3.0*H2+0.5*H3
  1051.       LU=SQRT(AA+CC)
  1052.       LV=SQRT(BB+DD)
  1053.       THXU=ATAN2(C,A)
  1054.       THUV=ATAN2(D,B)-THXU
  1055.       CSUV=COS(THUV)
  1056.       P41=5.0*LV*CSUV/LU*P50
  1057.       P14=5.0*LU*CSUV/LV*P05
  1058.       H1=ZV(2)-P01-P11-P41
  1059.       H2=ZUV(2)-P11-4.0*P41
  1060.       P21= 3.0*H1-H2
  1061.       P31=-2.0*H1+H2
  1062.       H1=ZU(3)-P10-P11-P14
  1063.       H2=ZUV(3)-P11-4.0*P14
  1064.       P12= 3.0*H1-H2
  1065.       P13=-2.0*H1+H2
  1066.       THUS=ATAN2(D-C,B-A)-THXU
  1067.       THSV=THUV-THUS
  1068.       AA= SIN(THSV)/LU
  1069.       BB=-COS(THSV)/LU
  1070.       CC= SIN(THUS)/LV
  1071.       DD= COS(THUS)/LV
  1072.       AC=AA*CC
  1073.       AD=AA*DD
  1074.       BC=BB*CC
  1075.       G1=AA*AC*(3.0*BC+2.0*AD)
  1076.       G2=CC*AC*(3.0*AD+2.0*BC)
  1077.       H1=-AA*AA*AA*(5.0*AA*BB*P50+(4.0*BC+AD)*P41)
  1078.      1   -CC*CC*CC*(5.0*CC*DD*P05+(4.0*AD+BC)*P14)
  1079.       H2=0.5*ZVV(2)-P02-P12
  1080.       H3=0.5*ZUU(3)-P20-P21
  1081.       P22=(G1*H2+G2*H3-H1)/(G1+G2)
  1082.       P32=H2-P22
  1083.       P23=H3-P22
  1084.       ITPV=IT0
  1085. C CONVERTS XII AND YII TO U-V SYSTEM.
  1086.    30 DX=XII-X0
  1087.       DY=YII-Y0
  1088.       U=AP*DX+BP*DY
  1089.       V=CP*DX+DP*DY
  1090. C EVALUATES THE POLYNOMIAL.
  1091.    31 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05))))
  1092.       P1=P10+V*(P11+V*(P12+V*(P13+V*P14)))
  1093.       P2=P20+V*(P21+V*(P22+V*P23))
  1094.       P3=P30+V*(P31+V*P32)
  1095.       P4=P40+V*P41
  1096.       ZII=P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5))))
  1097.       RETURN
  1098. C CALCULATION OF ZII BY EXTRAPOLATION IN THE RECTANGLE.
  1099. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
  1100.    40 IF(IT0.EQ.ITPV)     GO TO 50
  1101. C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE END
  1102. C POINTS OF THE BORDER LINE SEGMENT.
  1103.    41 JIPL=3*(IL1-1)
  1104.       JPD=0
  1105.       DO 43  I=1,2
  1106.         JIPL=JIPL+1
  1107.         IDP=IPL(JIPL)
  1108.         X(I)=XD(IDP)
  1109.         Y(I)=YD(IDP)
  1110.         Z(I)=ZD(IDP)
  1111.         JPDD=5*(IDP-1)
  1112.         DO 42  KPD=1,5
  1113.           JPD=JPD+1
  1114.           JPDD=JPDD+1
  1115.           PD(JPD)=PDD(JPDD)
  1116.    42   CONTINUE
  1117.    43 CONTINUE
  1118. C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
  1119. C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
  1120. C AND VICE VERSA.
  1121.    44 X0=X(1)
  1122.       Y0=Y(1)
  1123.       A=Y(2)-Y(1)
  1124.       B=X(2)-X(1)
  1125.       C=-B
  1126.       D=A
  1127.       AD=A*D
  1128.       BC=B*C
  1129.       DLT=AD-BC
  1130.       AP= D/DLT
  1131.       BP=-B/DLT
  1132.       CP=-BP
  1133.       DP= AP
  1134. C CONVERTS THE PARTIAL DERIVATIVES AT THE END POINTS OF THE
  1135. C BORDER LINE SEGMENT FOR THE U-V COORDINATE SYSTEM.
  1136.    45 AA=A*A
  1137.       ACT2=2.0*A*C
  1138.       CC=C*C
  1139.       AB=A*B
  1140.       ADBC=AD+BC
  1141.       CD=C*D
  1142.       BB=B*B
  1143.       BDT2=2.0*B*D
  1144.       DD=D*D
  1145.       DO 46  I=1,2
  1146.         JPD=5*I
  1147.         ZU(I)=A*PD(JPD-4)+C*PD(JPD-3)
  1148.         ZV(I)=B*PD(JPD-4)+D*PD(JPD-3)
  1149.         ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD)
  1150.         ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD)
  1151.         ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD)
  1152.    46 CONTINUE
  1153. C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
  1154.    47 P00=Z(1)
  1155.       P10=ZU(1)
  1156.       P01=ZV(1)
  1157.       P20=0.5*ZUU(1)
  1158.       P11=ZUV(1)
  1159.       P02=0.5*ZVV(1)
  1160.       H1=Z(2)-P00-P01-P02
  1161.       H2=ZV(2)-P01-ZVV(1)
  1162.       H3=ZVV(2)-ZVV(1)
  1163.       P03= 10.0*H1-4.0*H2+0.5*H3
  1164.       P04=-15.0*H1+7.0*H2    -H3
  1165.       P05=  6.0*H1-3.0*H2+0.5*H3
  1166.       H1=ZU(2)-P10-P11
  1167.       H2=ZUV(2)-P11
  1168.       P12= 3.0*H1-H2
  1169.       P13=-2.0*H1+H2
  1170.       P21=0.0
  1171.       P23=-ZUU(2)+ZUU(1)
  1172.       P22=-1.5*P23
  1173.       ITPV=IT0
  1174. C CONVERTS XII AND YII TO U-V SYSTEM.
  1175.    50 DX=XII-X0
  1176.       DY=YII-Y0
  1177.       U=AP*DX+BP*DY
  1178.       V=CP*DX+DP*DY
  1179. C EVALUATES THE POLYNOMIAL.
  1180.    51 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05))))
  1181.       P1=P10+V*(P11+V*(P12+V*P13))
  1182.       P2=P20+V*(P21+V*(P22+V*P23))
  1183.       ZII=P0+U*(P1+U*P2)
  1184.       RETURN
  1185. C CALCULATION OF ZII BY EXTRAPOLATION IN THE TRIANGLE.
  1186. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED.
  1187.    60 IF(IT0.EQ.ITPV)     GO TO 70
  1188. C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE VERTEX
  1189. C OF THE TRIANGLE.
  1190.    61 JIPL=3*IL2-2
  1191.       IDP=IPL(JIPL)
  1192.       X(1)=XD(IDP)
  1193.       Y(1)=YD(IDP)
  1194.       Z(1)=ZD(IDP)
  1195.       JPDD=5*(IDP-1)
  1196.       DO 62  KPD=1,5
  1197.         JPDD=JPDD+1
  1198.         PD(KPD)=PDD(JPDD)
  1199.    62 CONTINUE
  1200. C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL.
  1201.    67 P00=Z(1)
  1202.       P10=PD(1)
  1203.       P01=PD(2)
  1204.       P20=0.5*PD(3)
  1205.       P11=PD(4)
  1206.       P02=0.5*PD(5)
  1207.       ITPV=IT0
  1208. C CONVERTS XII AND YII TO U-V SYSTEM.
  1209.    70 U=XII-X(1)
  1210.       V=YII-Y(1)
  1211. C EVALUATES THE POLYNOMIAL.
  1212.    71 P0=P00+V*(P01+V*P02)
  1213.       P1=P10+V*P11
  1214.       ZII=P0+U*(P1+U*P20)
  1215.       RETURN
  1216.       END
  1217.       SUBROUTINE  IDSFFT(MD,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI,
  1218.      1                   IWK,WK)
  1219. C THIS SUBROUTINE PERFORMS SMOOTH SURFACE FITTING WHEN THE PRO-
  1220. C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
  1221. C DISTRIBUTED IN THE PLANE.
  1222. C THE INPUT PARAMETERS ARE
  1223. C     MD  = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
  1224. C         = 1 FOR NEW NCP AND/OR NEW XD-YD,
  1225. C         = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
  1226. C         = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
  1227. C     NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
  1228. C           MATING PARTIAL DERIVATIVES AT EACH DATA POINT
  1229. C           (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
  1230. C     NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
  1231. C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE X
  1232. C           COORDINATES OF THE DATA POINTS,
  1233. C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE Y
  1234. C           COORDINATES OF THE DATA POINTS,
  1235. C     ZD  = ARRAY OF DIMENSION NDP CONTAINING THE Z
  1236. C           COORDINATES OF THE DATA POINTS,
  1237. C     NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE
  1238. C           (MUST BE 1 OR GREATER),
  1239. C     NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE
  1240. C           (MUST BE 1 OR GREATER),
  1241. C     XI  = ARRAY OF DIMENSION NXI CONTAINING THE X
  1242. C           COORDINATES OF THE OUTPUT GRID POINTS,
  1243. C     YI  = ARRAY OF DIMENSION NYI CONTAINING THE Y
  1244. C           COORDINATES OF THE OUTPUT GRID POINTS.
  1245. C THE OUTPUT PARAMETER IS
  1246. C     ZI  = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (NXI,NYI),
  1247. C           WHERE THE INTERPOLATED Z VALUES AT THE OUTPUT
  1248. C           GRID POINTS ARE TO BE STORED.
  1249. C THE OTHER PARAMETERS ARE
  1250. C     IWK = INTEGER ARRAY OF DIMENSION
  1251. C              MAX0(31,27+NCP)*NDP+NXI*NYI
  1252. C           USED INTERNALLY AS A WORK AREA,
  1253. C     WK  = ARRAY OF DIMENSION 5*NDP USED INTERNALLY AS A
  1254. C           WORK AREA.
  1255. C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
  1256. C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
  1257. C YD ARRAYS MUST BE MADE WITH MD=1.  THE CALL WITH MD=2 MUST BE
  1258. C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
  1259. C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS.  THE CALL WITH
  1260. C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
  1261. C NXI, AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD, YD,
  1262. C XI, AND YI ARRAYS.  BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
  1263. C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
  1264. C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
  1265. C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
  1266. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
  1267. C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
  1268. C THEREFORE, SYSTEM DEPENDENT.
  1269. C THIS SUBROUTINE CALLS THE IDCLDP, IDGRID, IDPDRV, IDPTIP, AND
  1270. C IDTANG SUBROUTINES.
  1271. C DECLARATION STATEMENTS
  1272.       DIMENSION   XD(100),YD(100),ZD(100),XI(101),YI(101),
  1273.      1            ZI(10201),IWK(13301),WK(500)
  1274.       COMMON/IDPI/ITPV
  1275.       DATA  LUN/6/
  1276. C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
  1277. C (FOR MD=1,2,3)
  1278.    10 MD0=MD
  1279.       NCP0=NCP
  1280.       NDP0=NDP
  1281.       NXI0=NXI
  1282.       NYI0=NYI
  1283. C ERROR CHECK.  (FOR MD=1,2,3)
  1284.    20 IF(MD0.LT.1.OR.MD0.GT.3)           GO TO 90
  1285.       IF(NCP0.LT.2.OR.NCP0.GE.NDP0)      GO TO 90
  1286.       IF(NDP0.LT.4)                      GO TO 90
  1287.       IF(NXI0.LT.1.OR.NYI0.LT.1)         GO TO 90
  1288.       IF(MD0.GE.2)        GO TO 21
  1289.       IWK(1)=NCP0
  1290.       IWK(2)=NDP0
  1291.       GO TO 22
  1292.    21 NCPPV=IWK(1)
  1293.       NDPPV=IWK(2)
  1294.       IF(NCP0.NE.NCPPV)   GO TO 90
  1295.       IF(NDP0.NE.NDPPV)   GO TO 90
  1296.    22 IF(MD0.GE.3)        GO TO 23
  1297.       IWK(3)=NXI0
  1298.       IWK(4)=NYI0
  1299.       GO TO 30
  1300.    23 NXIPV=IWK(3)
  1301.       NYIPV=IWK(4)
  1302.       IF(NXI0.NE.NXIPV)   GO TO 90
  1303.       IF(NYI0.NE.NYIPV)   GO TO 90
  1304. C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY.  (FOR MD=1,2,3)
  1305.    30 JWIPT=16
  1306.       JWIWL=6*NDP0+1
  1307.       JWNGP0=JWIWL-1
  1308.       JWIPL=24*NDP0+1
  1309.       JWIWP=30*NDP0+1
  1310.       JWIPC=27*NDP0+1
  1311.       JWIGP0=MAX0(31,27+NCP0)*NDP0
  1312. C TRIANGULATES THE X-Y PLANE.  (FOR MD=1)
  1313.    40 IF(MD0.GT.1)   GO TO 50
  1314.       CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
  1315.      1            IWK(JWIWL),IWK(JWIWP),WK)
  1316.       IWK(5)=NT
  1317.       IWK(6)=NL
  1318.       IF(NT.EQ.0)    RETURN
  1319. C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT.  (FOR MD=1)
  1320.    50 IF(MD0.GT.1)   GO TO 60
  1321.       CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC))
  1322.       IF(IWK(JWIPC).EQ.0)      RETURN
  1323. C SORTS OUTPUT GRID POINTS IN ASCENDING ORDER OF THE TRIANGLE
  1324. C NUMBER AND THE BORDER LINE SEGMENT NUMBER.  (FOR MD=1,2)
  1325.    60 IF(MD0.EQ.3)   GO TO 70
  1326.       CALL IDGRID(XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),NXI0,NYI0,
  1327.      1            XI,YI,IWK(JWNGP0+1),IWK(JWIGP0+1))
  1328. C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
  1329. C (FOR MD=1,2,3)
  1330.    70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
  1331. C INTERPOLATES THE ZI VALUES.  (FOR MD=1,2,3)
  1332.    80 ITPV=0
  1333.       JIG0MX=0
  1334.       JIG1MN=NXI0*NYI0+1
  1335.       NNGP=NT+2*NL
  1336.       DO 89  JNGP=1,NNGP
  1337.         ITI=JNGP
  1338.         IF(JNGP.LE.NT)    GO TO 81
  1339.         IL1=(JNGP-NT+1)/2
  1340.         IL2=(JNGP-NT+2)/2
  1341.         IF(IL2.GT.NL)     IL2=1
  1342.         ITI=IL1*(NT+NL)+IL2
  1343.    81   JWNGP=JWNGP0+JNGP
  1344.         NGP0=IWK(JWNGP)
  1345.         IF(NGP0.EQ.0)     GO TO 86
  1346.         JIG0MN=JIG0MX+1
  1347.         JIG0MX=JIG0MX+NGP0
  1348.         DO 82  JIGP=JIG0MN,JIG0MX
  1349.           JWIGP=JWIGP0+JIGP
  1350.           IZI=IWK(JWIGP)
  1351.           IYI=(IZI-1)/NXI0+1
  1352.           IXI=IZI-NXI0*(IYI-1)
  1353.           CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
  1354.      1                ITI,XI(IXI),YI(IYI),ZI(IZI))
  1355.    82   CONTINUE
  1356.    86   JWNGP=JWNGP0+2*NNGP+1-JNGP
  1357.         NGP1=IWK(JWNGP)
  1358.         IF(NGP1.EQ.0)     GO TO 89
  1359.         JIG1MX=JIG1MN-1
  1360.         JIG1MN=JIG1MN-NGP1
  1361.         DO 87  JIGP=JIG1MN,JIG1MX
  1362.           JWIGP=JWIGP0+JIGP
  1363.           IZI=IWK(JWIGP)
  1364.           IYI=(IZI-1)/NXI0+1
  1365.           IXI=IZI-NXI0*(IYI-1)
  1366.           CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
  1367.      1                ITI,XI(IXI),YI(IYI),ZI(IZI))
  1368.    87   CONTINUE
  1369.    89 CONTINUE
  1370.       RETURN
  1371. C ERROR EXIT
  1372.    90 WRITE (LUN,2090) MD0,NCP0,NDP0,NXI0,NYI0
  1373.       RETURN
  1374. C FORMAT STATEMENT FOR ERROR MESSAGE
  1375.  2090 FORMAT(1X/41H ***   IMPROPER INPUT PARAMETER VALUE(S)./
  1376.      1   7H   MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6,
  1377.      2   10X,5HNXI =,I6,10X,5HNYI =,I6/
  1378.      3   35H ERROR DETECTED IN ROUTINE   IDSFFT/)
  1379.       END
  1380.       SUBROUTINE  IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
  1381. C THIS SUBROUTINE PERFORMS TRIANGULATION.  IT DIVIDES THE X-Y
  1382. C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA
  1383. C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE
  1384. C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS
  1385. C CORRESPONDING TO THE BORDER LINE SEGMENTS.
  1386. C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE
  1387. C ARE LISTED COUNTER-CLOCKWISE.  POINT NUMBERS OF THE END POINTS
  1388. C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE,
  1389. C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.
  1390. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
  1391. C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
  1392. C THEREFORE, SYSTEM DEPENDENT.
  1393. C THIS SUBROUTINE CALLS THE IDXCHG FUNCTION.
  1394. C THE INPUT PARAMETERS ARE
  1395. C     NDP = NUMBER OF DATA POINTS,
  1396. C     XD  = ARRAY OF DIMENSION NDP CONTAINING THE
  1397. C           X COORDINATES OF THE DATA POINTS,
  1398. C     YD  = ARRAY OF DIMENSION NDP CONTAINING THE
  1399. C           Y COORDINATES OF THE DATA POINTS.
  1400. C THE OUTPUT PARAMETERS ARE
  1401. C     NT  = NUMBER OF TRIANGLES,
  1402. C     IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE
  1403. C           POINT NUMBERS OF THE VERTEXES OF THE (IT)TH
  1404. C           TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND,
  1405. C           (3*IT-1)ST, AND (3*IT)TH ELEMENTS,
  1406. C           IT=1,2,...,NT,
  1407. C     NL  = NUMBER OF BORDER LINE SEGMENTS,
  1408. C     IPL = INTEGER ARRAY OF DIMENSION 6*NDP, WHERE THE
  1409. C           POINT NUMBERS OF THE END POINTS OF THE (IL)TH
  1410. C           BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE
  1411. C           NUMBER ARE TO BE STORED AS THE (3*IL-2)ND,
  1412. C           (3*IL-1)ST, AND (3*IL)TH ELEMENTS,
  1413. C           IL=1,2,..., NL.
  1414. C THE OTHER PARAMETERS ARE
  1415. C     IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
  1416. C           INTERNALLY AS A WORK AREA,
  1417. C     IWP = INTEGER ARRAY OF DIMENSION NDP USED
  1418. C           INTERNALLY AS A WORK AREA,
  1419. C     WK  = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
  1420. C           WORK AREA.
  1421. C DECLARATION STATEMENTS
  1422.       DIMENSION   XD(100),YD(100),IPT(585),IPL(600),
  1423.      1            IWL(1800),IWP(100),WK(100)
  1424.       DIMENSION   ITF(2)
  1425.       DATA  RATIO/1.0E-6/, NREP/100/, LUN/6/
  1426. C STATEMENT FUNCTIONS
  1427.       DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
  1428.       SIDE(U1,V1,U2,V2,U3,V3)=(V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
  1429. C PRELIMINARY PROCESSING
  1430.    10 NDP0=NDP
  1431.       NDPM1=NDP0-1
  1432.       IF(NDP0.LT.4)       GO TO 90
  1433. C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
  1434.    20 DSQMN=DSQF(XD(1),YD(1),XD(2),YD(2))
  1435.       IPMN1=1
  1436.       IPMN2=2
  1437.       DO 22  IP1=1,NDPM1
  1438.         X1=XD(IP1)
  1439.         Y1=YD(IP1)
  1440.         IP1P1=IP1+1
  1441.         DO 21  IP2=IP1P1,NDP0
  1442.           DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
  1443.           IF(DSQI.EQ.0.0)      GO TO 91
  1444.           IF(DSQI.GE.DSQMN)    GO TO 21
  1445.           DSQMN=DSQI
  1446.           IPMN1=IP1
  1447.           IPMN2=IP2
  1448.    21   CONTINUE
  1449.    22 CONTINUE
  1450.       DSQ12=DSQMN
  1451.       XDMP=(XD(IPMN1)+XD(IPMN2))/2.0
  1452.       YDMP=(YD(IPMN1)+YD(IPMN2))/2.0
  1453. C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
  1454. C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
  1455. C NUMBERS IN THE IWP ARRAY.
  1456.    30 JP1=2
  1457.       DO 31  IP1=1,NDP0
  1458.         IF(IP1.EQ.IPMN1.OR.IP1.EQ.IPMN2)      GO TO 31
  1459.         JP1=JP1+1
  1460.         IWP(JP1)=IP1
  1461.         WK(JP1)=DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
  1462.    31 CONTINUE
  1463.       DO 33  JP1=3,NDPM1
  1464.         DSQMN=WK(JP1)
  1465.         JPMN=JP1
  1466.         DO 32  JP2=JP1,NDP0
  1467.           IF(WK(JP2).GE.DSQMN)      GO TO 32
  1468.           DSQMN=WK(JP2)
  1469.           JPMN=JP2
  1470.    32   CONTINUE
  1471.         ITS=IWP(JP1)
  1472.         IWP(JP1)=IWP(JPMN)
  1473.         IWP(JPMN)=ITS
  1474.         WK(JPMN)=WK(JP1)
  1475.    33 CONTINUE
  1476. C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
  1477. C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
  1478.    35 AR=DSQ12*RATIO
  1479.       X1=XD(IPMN1)
  1480.       Y1=YD(IPMN1)
  1481.       DX21=XD(IPMN2)-X1
  1482.       DY21=YD(IPMN2)-Y1
  1483.       DO 36  JP=3,NDP0
  1484.         IP=IWP(JP)
  1485.         IF(ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21).GT.AR)
  1486.      1               GO TO 37
  1487.    36 CONTINUE
  1488.       GO TO 92
  1489.    37 IF(JP.EQ.3)    GO TO 40
  1490.       JPMX=JP
  1491.       JP=JPMX+1
  1492.       DO 38  JPC=4,JPMX
  1493.         JP=JP-1
  1494.         IWP(JP)=IWP(JP-1)
  1495.    38 CONTINUE
  1496.       IWP(3)=IP
  1497. C FORMS THE FIRST TRIANGLE.  STORES POINT NUMBERS OF THE VER-
  1498. C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
  1499. C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
  1500. C THE IPL ARRAY.
  1501.    40 IP1=IPMN1
  1502.       IP2=IPMN2
  1503.       IP3=IWP(3)
  1504.       IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
  1505.      1     .GE.0.0)       GO TO 41
  1506.       IP1=IPMN2
  1507.       IP2=IPMN1
  1508.    41 NT0=1
  1509.       NTT3=3
  1510.       IPT(1)=IP1
  1511.       IPT(2)=IP2
  1512.       IPT(3)=IP3
  1513.       NL0=3
  1514.       NLT3=9
  1515.       IPL(1)=IP1
  1516.       IPL(2)=IP2
  1517.       IPL(3)=1
  1518.       IPL(4)=IP2
  1519.       IPL(5)=IP3
  1520.       IPL(6)=1
  1521.       IPL(7)=IP3
  1522.       IPL(8)=IP1
  1523.       IPL(9)=1
  1524. C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
  1525.    50 DO 79  JP1=4,NDP0
  1526.         IP1=IWP(JP1)
  1527.         X1=XD(IP1)
  1528.         Y1=YD(IP1)
  1529. C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
  1530.         IP2=IPL(1)
  1531.         JPMN=1
  1532.         DXMN=XD(IP2)-X1
  1533.         DYMN=YD(IP2)-Y1
  1534.         DSQMN=DXMN**2+DYMN**2
  1535.         ARMN=DSQMN*RATIO
  1536.         JPMX=1
  1537.         DXMX=DXMN
  1538.         DYMX=DYMN
  1539.         DSQMX=DSQMN
  1540.         ARMX=ARMN
  1541.         DO 52  JP2=2,NL0
  1542.           IP2=IPL(3*JP2-2)
  1543.           DX=XD(IP2)-X1
  1544.           DY=YD(IP2)-Y1
  1545.           AR=DY*DXMN-DX*DYMN
  1546.           IF(AR.GT.ARMN)       GO TO 51
  1547.           DSQI=DX**2+DY**2
  1548.           IF(AR.GE.(-ARMN).AND.DSQI.GE.DSQMN)      GO TO 51
  1549.           JPMN=JP2
  1550.           DXMN=DX
  1551.           DYMN=DY
  1552.           DSQMN=DSQI
  1553.           ARMN=DSQMN*RATIO
  1554.    51     AR=DY*DXMX-DX*DYMX
  1555.           IF(AR.LT.(-ARMX))    GO TO 52
  1556.           DSQI=DX**2+DY**2
  1557.           IF(AR.LE.ARMX.AND.DSQI.GE.DSQMX)    GO TO 52
  1558.           JPMX=JP2
  1559.           DXMX=DX
  1560.           DYMX=DY
  1561.           DSQMX=DSQI
  1562.           ARMX=DSQMX*RATIO
  1563.    52   CONTINUE
  1564.         IF(JPMX.LT.JPMN)  JPMX=JPMX+NL0
  1565.         NSH=JPMN-1
  1566.         IF(NSH.LE.0)      GO TO 60
  1567. C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
  1568. C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
  1569.         NSHT3=NSH*3
  1570.         DO 53  JP2T3=3,NSHT3,3
  1571.           JP3T3=JP2T3+NLT3
  1572.           IPL(JP3T3-2)=IPL(JP2T3-2)
  1573.           IPL(JP3T3-1)=IPL(JP2T3-1)
  1574.           IPL(JP3T3)  =IPL(JP2T3)
  1575.    53   CONTINUE
  1576.         DO 54  JP2T3=3,NLT3,3
  1577.           JP3T3=JP2T3+NSHT3
  1578.           IPL(JP2T3-2)=IPL(JP3T3-2)
  1579.           IPL(JP2T3-1)=IPL(JP3T3-1)
  1580.           IPL(JP2T3)  =IPL(JP3T3)
  1581.    54   CONTINUE
  1582.         JPMX=JPMX-NSH
  1583. C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
  1584. C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
  1585. C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
  1586.    60   JWL=0
  1587.         DO 64  JP2=JPMX,NL0
  1588.           JP2T3=JP2*3
  1589.           IPL1=IPL(JP2T3-2)
  1590.           IPL2=IPL(JP2T3-1)
  1591.           IT  =IPL(JP2T3)
  1592. C - - ADDS A TRIANGLE TO THE IPT ARRAY.
  1593.           NT0=NT0+1
  1594.           NTT3=NTT3+3
  1595.           IPT(NTT3-2)=IPL2
  1596.           IPT(NTT3-1)=IPL1
  1597.           IPT(NTT3)  =IP1
  1598. C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
  1599.           IF(JP2.NE.JPMX)      GO TO 61
  1600.           IPL(JP2T3-1)=IP1
  1601.           IPL(JP2T3)  =NT0
  1602.    61     IF(JP2.NE.NL0)       GO TO 62
  1603.           NLN=JPMX+1
  1604.           NLNT3=NLN*3
  1605.           IPL(NLNT3-2)=IP1
  1606.           IPL(NLNT3-1)=IPL(1)
  1607.           IPL(NLNT3)  =NT0
  1608. C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
  1609. C - - LINE SEGMENTS.
  1610.    62     ITT3=IT*3
  1611.           IPTI=IPT(ITT3-2)
  1612.           IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2)   GO TO 63
  1613.           IPTI=IPT(ITT3-1)
  1614.           IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2)   GO TO 63
  1615.           IPTI=IPT(ITT3)
  1616. C - - CHECKS IF THE EXCHANGE IS NECESSARY.
  1617.    63     IF(IDXCHG(XD,YD,IP1,IPTI,IPL1,IPL2).EQ.0)     GO TO 64
  1618. C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
  1619.           IPT(ITT3-2)=IPTI
  1620.           IPT(ITT3-1)=IPL1
  1621.           IPT(ITT3)  =IP1
  1622.           IPT(NTT3-1)=IPTI
  1623.           IF(JP2.EQ.JPMX)      IPL(JP2T3)=IT
  1624.           IF(JP2.EQ.NL0.AND.IPL(3).EQ.IT)     IPL(3)=NT0
  1625. C - - SETS FLAGS IN THE IWL ARRAY.
  1626.           JWL=JWL+4
  1627.           IWL(JWL-3)=IPL1
  1628.           IWL(JWL-2)=IPTI
  1629.           IWL(JWL-1)=IPTI
  1630.           IWL(JWL)  =IPL2
  1631.    64   CONTINUE
  1632.         NL0=NLN
  1633.         NLT3=NLNT3
  1634.         NLF=JWL/2
  1635.         IF(NLF.EQ.0)      GO TO 79
  1636. C - IMPROVES TRIANGULATION.
  1637.    70   NTT3P3=NTT3+3
  1638.         DO 78  IREP=1,NREP
  1639.           DO 76  ILF=1,NLF
  1640.             ILFT2=ILF*2
  1641.             IPL1=IWL(ILFT2-1)
  1642.             IPL2=IWL(ILFT2)
  1643. C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
  1644. C - - THE FLAGGED LINE SEGMENT.
  1645.             NTF=0
  1646.             DO 71  ITT3R=3,NTT3,3
  1647.               ITT3=NTT3P3-ITT3R
  1648.               IPT1=IPT(ITT3-2)
  1649.               IPT2=IPT(ITT3-1)
  1650.               IPT3=IPT(ITT3)
  1651.               IF(IPL1.NE.IPT1.AND.IPL1.NE.IPT2.AND.
  1652.      1           IPL1.NE.IPT3)      GO TO 71
  1653.               IF(IPL2.NE.IPT1.AND.IPL2.NE.IPT2.AND.
  1654.      1           IPL2.NE.IPT3)      GO TO 71
  1655.               NTF=NTF+1
  1656.               ITF(NTF)=ITT3/3
  1657.               IF(NTF.EQ.2)     GO TO 72
  1658.    71       CONTINUE
  1659.             IF(NTF.LT.2)       GO TO 76
  1660. C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
  1661. C - - ON THE LINE SEGMENT.
  1662.    72       IT1T3=ITF(1)*3
  1663.             IPTI1=IPT(IT1T3-2)
  1664.             IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2)    GO TO 73
  1665.             IPTI1=IPT(IT1T3-1)
  1666.             IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2)    GO TO 73
  1667.             IPTI1=IPT(IT1T3)
  1668.    73       IT2T3=ITF(2)*3
  1669.             IPTI2=IPT(IT2T3-2)
  1670.             IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2)    GO TO 74
  1671.             IPTI2=IPT(IT2T3-1)
  1672.             IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2)    GO TO 74
  1673.             IPTI2=IPT(IT2T3)
  1674. C - - CHECKS IF THE EXCHANGE IS NECESSARY.
  1675.    74       IF(IDXCHG(XD,YD,IPTI1,IPTI2,IPL1,IPL2).EQ.0)
  1676.      1         GO TO 76
  1677. C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
  1678.             IPT(IT1T3-2)=IPTI1
  1679.             IPT(IT1T3-1)=IPTI2
  1680.             IPT(IT1T3)  =IPL1
  1681.             IPT(IT2T3-2)=IPTI2
  1682.             IPT(IT2T3-1)=IPTI1
  1683.             IPT(IT2T3)  =IPL2
  1684. C - - SETS NEW FLAGS.
  1685.             JWL=JWL+8
  1686.             IWL(JWL-7)=IPL1
  1687.             IWL(JWL-6)=IPTI1
  1688.             IWL(JWL-5)=IPTI1
  1689.             IWL(JWL-4)=IPL2
  1690.             IWL(JWL-3)=IPL2
  1691.             IWL(JWL-2)=IPTI2
  1692.             IWL(JWL-1)=IPTI2
  1693.             IWL(JWL)  =IPL1
  1694.             DO 75  JLT3=3,NLT3,3
  1695.               IPLJ1=IPL(JLT3-2)
  1696.               IPLJ2=IPL(JLT3-1)
  1697.               IF((IPLJ1.EQ.IPL1.AND.IPLJ2.EQ.IPTI2).OR.
  1698.      1           (IPLJ2.EQ.IPL1.AND.IPLJ1.EQ.IPTI2))
  1699.      2                         IPL(JLT3)=ITF(1)
  1700.               IF((IPLJ1.EQ.IPL2.AND.IPLJ2.EQ.IPTI1).OR.
  1701.      1           (IPLJ2.EQ.IPL2.AND.IPLJ1.EQ.IPTI1))
  1702.      2                         IPL(JLT3)=ITF(2)
  1703.    75       CONTINUE
  1704.    76     CONTINUE
  1705.           NLFC=NLF
  1706.           NLF=JWL/2
  1707.           IF(NLF.EQ.NLFC)      GO TO 79
  1708. C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
  1709.           JWL=0
  1710.           JWL1MN=(NLFC+1)*2
  1711.           NLFT2=NLF*2
  1712.           DO 77  JWL1=JWL1MN,NLFT2,2
  1713.             JWL=JWL+2
  1714.             IWL(JWL-1)=IWL(JWL1-1)
  1715.             IWL(JWL)  =IWL(JWL1)
  1716.    77     CONTINUE
  1717.           NLF=JWL/2
  1718.    78   CONTINUE
  1719.    79 CONTINUE
  1720. C REARRANGES THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
  1721. C ARE LISTED COUNTER-CLOCKWISE.
  1722.    80 DO 81  ITT3=3,NTT3,3
  1723.         IP1=IPT(ITT3-2)
  1724.         IP2=IPT(ITT3-1)
  1725.         IP3=IPT(ITT3)
  1726.         IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
  1727.      1       .GE.0.0)     GO TO 81
  1728.         IPT(ITT3-2)=IP2
  1729.         IPT(ITT3-1)=IP1
  1730.    81 CONTINUE
  1731.       NT=NT0
  1732.       NL=NL0
  1733.       RETURN
  1734. C ERROR EXIT
  1735.    90 WRITE (LUN,2090)  NDP0
  1736.       GO TO 93
  1737.    91 WRITE (LUN,2091)  NDP0,IP1,IP2,X1,Y1
  1738.       GO TO 93
  1739.    92 WRITE (LUN,2092)  NDP0
  1740.    93 WRITE (LUN,2093)
  1741.       NT=0
  1742.       RETURN
  1743. C FORMAT STATEMENTS
  1744.  2090 FORMAT(1X/23H ***   NDP LESS THAN 4./8H   NDP =,I5)
  1745.  2091 FORMAT(1X/29H ***   IDENTICAL DATA POINTS./
  1746.      1   8H   NDP =,I5,5X,5HIP1 =,I5,5X,5HIP2 =,I5,
  1747.      2   5X,4HXD =,E12.4,5X,4HYD =,E12.4)
  1748.  2092 FORMAT(1X/33H ***   ALL COLLINEAR DATA POINTS./
  1749.      1   8H   NDP =,I5)
  1750.  2093 FORMAT(35H ERROR DETECTED IN ROUTINE   IDTANG/)
  1751.       END
  1752.       FUNCTION  IDXCHG(X,Y,I1,I2,I3,I4)
  1753. C THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO
  1754. C TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION
  1755. C BY C. L. LAWSON.
  1756. C THE INPUT PARAMETERS ARE
  1757. C     X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA
  1758. C           POINTS,
  1759. C     I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1, P2,
  1760. C           P3, AND P4 THAT FORM A QUADRILATERAL WITH P3
  1761. C           AND P4 CONNECTED DIAGONALLY.
  1762. C THIS FUNCTION RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EX-
  1763. C CHANGE IS NECESSARY, AND 0 (ZERO) OTHERWISE.
  1764. C DECLARATION STATEMENTS
  1765.       DIMENSION   X(100),Y(100)
  1766.       EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ),
  1767.      1            (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ)
  1768. C PRELIMINARY PROCESSING
  1769.    10 X1=X(I1)
  1770.       Y1=Y(I1)
  1771.       X2=X(I2)
  1772.       Y2=Y(I2)
  1773.       X3=X(I3)
  1774.       Y3=Y(I3)
  1775.       X4=X(I4)
  1776.       Y4=Y(I4)
  1777. C CALCULATION
  1778.    20 IDX=0
  1779.       U3=(Y2-Y3)*(X1-X3)-(X2-X3)*(Y1-Y3)
  1780.       U4=(Y1-Y4)*(X2-X4)-(X1-X4)*(Y2-Y4)
  1781.       IF(U3*U4.LE.0.0)    GO TO 30
  1782.       U1=(Y3-Y1)*(X4-X1)-(X3-X1)*(Y4-Y1)
  1783.       U2=(Y4-Y2)*(X3-X2)-(X4-X2)*(Y3-Y2)
  1784.       A1SQ=(X1-X3)**2+(Y1-Y3)**2
  1785.       B1SQ=(X4-X1)**2+(Y4-Y1)**2
  1786.       C1SQ=(X3-X4)**2+(Y3-Y4)**2
  1787.       A2SQ=(X2-X4)**2+(Y2-Y4)**2
  1788.       B2SQ=(X3-X2)**2+(Y3-Y2)**2
  1789.       C3SQ=(X2-X1)**2+(Y2-Y1)**2
  1790.       S1SQ=U1*U1/(C1SQ*AMAX1(A1SQ,B1SQ))
  1791.       S2SQ=U2*U2/(C2SQ*AMAX1(A2SQ,B2SQ))
  1792.       S3SQ=U3*U3/(C3SQ*AMAX1(A3SQ,B3SQ))
  1793.       S4SQ=U4*U4/(C4SQ*AMAX1(A4SQ,B4SQ))
  1794.       IF(AMIN1(S1SQ,S2SQ).LT.AMIN1(S3SQ,S4SQ))     IDX=1
  1795.    30 IDXCHG=IDX
  1796.       RETURN
  1797.       END
  1798.  
  1799.